#---- function ----------
linear.detrend=function(temp){
	xx=1:length(temp)
	MAT=na.omit(cbind(xx,temp))
	spdata=array(NA,c(length(xx)))
	if (length(MAT)>=8){
		fit<-lm(MAT[,2]~MAT[,1])
		yy=xx*fit$coefficients[2]+fit$coefficients[1]
		spdata=temp-yy
	}
	return(spdata)
}

mov.detrend=function(temp){
	  ap=mov.avg.wgt(temp, c(1,2,3,2,1))
	  ap=mov.avg.wgt(temp, c(1,1,1,1,1))	  
	  return(as.numeric(temp-ap))
}

find.cor2=function(spdata,index,p_value){
	new=array(NA,dim(spdata))
	for (i in 1:dim(new)[1]){
		for (j in 1:dim(new)[2]){
			if (sum(!is.na(spdata[i,j,]))>=10){
			new[i,j,]=index
			}
		}
	}
	corr=find.cor(spdata,new,plim=p_value)
	return(corr)
}

test_siginificance=function(spdata,mu,plim){
significance=array(0,c(dim(spdata)[1],dim(spdata)[2]))
for (i in 1:dim(spdata)[1]){
	for (j in 1:dim(spdata)[2]){
		y=spdata[i,j,]
		if (sum(!is.na(y))>5){
			ap=t.test(y,mu=mu, alternative="two.sided")$p.value
			if (ap<plim) significance[i,j]=1
		}
	}	
}
return(significance)	
}

#---- read the gph5 data ---------
read_US=function(name,pressure,start_month,end_month){
date=array(NA,c(792,2))
for (i in 1:792){
  if (i%%12>0) {date[i,1]=1948+floor(i/12);date[i,2]=i%%12}
  if (i%%12==0){date[i,1]=1948+floor(i/12)-1;date[i,2]=12}
}	
datadir='Data/NCEP/'
datafile = open.ncdf(paste(datadir,name,sep=''))
lon = get.var.ncdf(datafile, varid='lon')
lat = get.var.ncdf(datafile, varid='lat')
ind1=which(lon>=220 & lon<=320)
ind2=which(lat>=20 & lat<=60)
ind4=which(date[,1]>=1999 & date[,1]<=2013)
lon=lon[ind1]
lat=lat[ind2]
date=date[ind4,]
level = get.var.ncdf(datafile,varid='level')
ind3=which(level== pressure)
origin_data = get.var.ncdf(datafile,start=c(ind1[1],ind2[1],ind3,ind4[1]), count=c(length(ind1),length(ind2),1,length(ind4)))
origin =origin_data[,length(lat):1,]
lat=rev(lat)
close.ncdf(datafile)
temp=array(NA,c(length(lon),length(lat),15))
for (year in 1999:2013){
  ind=(date[,1]==year & date[,2]>=start_month & date[,2]<=end_month)
  temp[,,year-1998]=apply(origin[,,ind],c(1,2),mean,na.rm=TRUE)
}
for (i in 1:dim(temp)[1]){
	for (j in 1:dim(temp)[2]){
		temp[i,j,]= mov.detrend(temp[i,j,])
	}
}
return(list('data'=temp,'lon'=lon,'lat'=lat))
}

#------ read suface ---------
read_surface=function(name,start_month,end_month){
date=array(NA,c(792,2))
for (i in 1:792){
  if (i%%12>0) {date[i,1]=1948+floor(i/12);date[i,2]=i%%12}
  if (i%%12==0){date[i,1]=1948+floor(i/12)-1;date[i,2]=12}
}
datadir='Data/NCEP/'	
datafile = open.ncdf(paste(datadir,name,sep=''))
lon = get.var.ncdf(datafile, varid='lon')
lat = get.var.ncdf(datafile, varid='lat')
ind1=which(lon>=220 & lon<=320)
ind2=which(lat>=20 & lat<=60)
ind4=which(date[,1]>=1999 & date[,1]<=2013)
lon=lon[ind1]
lat=lat[ind2]
date=date[ind4,]
origin_data = get.var.ncdf(datafile,start=c(ind1[1],ind2[1],ind4[1]), count=c(length(ind1),length(ind2),length(ind4)))
origin =origin_data[,length(lat):1,]
lat=rev(lat)
close.ncdf(datafile)

temp=array(NA,c(length(lon),length(lat),15))
for (year in 1999:2013){
  ind=(date[,1]==year & date[,2]>=start_month & date[,2]<=end_month)
  temp[,,year-1998]=apply(origin[,,ind],c(1,2),mean,na.rm=TRUE)
}
for (i in 1:dim(temp)[1]){
	for (j in 1:dim(temp)[2]){
		temp[i,j,]= mov.detrend(temp[i,j,])
	}
}
return(list('data'=temp,'lon'=lon,'lat'=lat))
}


read_met=function(name){
	date=array(NA,c(792,2))
	for (i in 1:792){
	  if (i%%12>0) {date[i,1]=1948+floor(i/12);date[i,2]=i%%12}
	  if (i%%12==0){date[i,1]=1948+floor(i/12)-1;date[i,2]=12}
	}
	datadir='Data/NCEP/'	
	datafile = open.ncdf(paste(datadir,name,sep=''))
	lon = get.var.ncdf(datafile, varid='lon')
	lat = get.var.ncdf(datafile, varid='lat')
	ind1=which(lon>=220 & lon<=320)
	ind2=which(lat>=20 & lat<=60)
	ind4=which(date[,1]>=1999 & date[,1]<=2013)
	lon=lon[ind1]
	lat=lat[ind2]
	date=date[ind4,]
	origin_data = get.var.ncdf(datafile,start=c(ind1[1],ind2[1],ind4[1]), 	count=c(length(ind1),length(ind2),length(ind4)))
	origin =origin_data[,length(lat):1,]
	lat=rev(lat)
	close.ncdf(datafile)
	
	spdata=array(NA,dim(origin))
	for (month in 1:12){
		  ind=(date[,2]==month)
		  y=origin[,,ind]
		  for (i in 1:dim(y)[1]){
			  for (j in 1:dim(y)[2]){
				  spdata[i,j,ind]= mov.detrend(y[i,j,])
			  }
		  }
	  }	
	return(list('data'= spdata,'lon'=lon,'lat'=lat))
}